    clear;
    close all
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % main table with welfare comparisons 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % degenerate distributions -- panel a)-g)
    load('./results/state1.mat');
    table_main = simulation_output.maintable;
    for i=2:7
     file = sprintf('./results/state%d.mat',i);
     load(file);
     table_main = cat(1,table_main,simulation_output.maintable);     
     clearvars simulation_output
    end

    % add panel h -- main results
    load('./results/allstates_main.mat')
    table_main = cat(1,table_main,simulation_output.maintable);  
    % keep only reported columns
    table_main.gainsGHHW=[];
    table_main.gainsGLTHI=[];
    table_main.metric1= [];
    writetable(table_main,'../../figsandtabs/table_main.xlsx');
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CALIBRATED PREMIUMS
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    clear;
    
    load('./results/allstates_main.mat')
    load('./results/model_parameters_main.mat')
    MU = model_parameters.MU;
    T = size(MU,3);
    age = [25:(25+T-1)]';
    
    Pspot_N = simulation_output.Pspot_N;
    Ppaid_N = simulation_output.Ppaid_N;
    Wio_abi= simulation_output.Wio_abi;
    Wio_rea= simulation_output.Wio_rea;
    Cg_abi = simulation_output.Cg_abi;
    Cg_rea = simulation_output.Cg_rea;
    iX = simulation_output.iX;
    Nstates = size(Pspot_N,1)-1;
    Ppaid_spot = simulation_output.Ppaid_spot;
    nInd = model_parameters.nInd;

    writetable(array2table(Pspot_N),'./results/Pspot_N_statexT.xlsx');
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%% CONSUMPTION %%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
      
    %%%%%%%%%%%%%%%%
    % Ed 13 (abi)
    %%%%%%%%%%%%%%%%
      
      % ST
      
      Cons_ST_abi =Wio_abi-Ppaid_spot;
      Cons_ST_abi_mean = sum(Cons_ST_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
      Cons_ST_abi_sd=sqrt(sum(((Cons_ST_abi-repmat(Cons_ST_abi_mean,1,nInd)).*(iX~=Nstates+1)).^2,2)./sum(iX~=(Nstates+1),2));

      % HHW
      
      Cons_HHW_abi =Consumption_HHW(Cg_abi,iX);
      Cons_HHW_abi_mean = sum(Cons_HHW_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
  
      % PHI
      
      Cons_PHI_abi =Wio_abi-Ppaid_N;
      Cons_PHI_abi_mean = sum(Cons_PHI_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
        
      f = gcf;
      plot(age,Cons_HHW_abi_mean,age,Cons_PHI_abi_mean,age,Cons_ST_abi_mean)
            
      xlabel('age', 'FontSize',15)
      ylabel('Avg. Consumption (Th. USD)','Fontsize',15)
      h = legend('HHW','PHI','ST');
      set(h,'Location','northwest');
      % b,r,k,g
      set(h,'fontsize',14);
      set(gca,'fontsize',14);

      F_hhw_phi_ST_abi_mean = [age,Cons_HHW_abi_mean,Cons_PHI_abi_mean,Cons_ST_abi_mean];
      writetable(array2table(F_hhw_phi_ST_abi_mean),'./results/F_hhw_phi_ST_abi_mean_HHW_PHI_ST.xlsx')    
             
      %%%%%%%%%%%% 
      % Ed 10 (rea)   
      %%%%%%%%%%%%% 
      
      % ST
      
      Cons_ST_rea =Wio_rea-Ppaid_spot;
      Cons_ST_rea_mean = sum(Cons_ST_rea.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 

      % HHW
      
      Cons_HHW_rea =Consumption_HHW(Cg_rea,iX);
      Cons_HHW_rea_mean = sum(Cons_HHW_rea.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
  
      % PHI
      
      Cons_PHI_rea =Wio_rea-Ppaid_N;
      Cons_PHI_rea_mean = sum(Cons_PHI_rea.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
     
      plot(age,Cons_HHW_rea_mean,age,Cons_PHI_rea_mean,age,Cons_ST_rea_mean)
     
      
      xlabel('age', 'FontSize',15)
      ylabel('Avg. Consumption (Th. USD)','Fontsize',15)
      h = legend('HHW','PHI','ST');
      set(h,'Location','northwest');
      set(h,'fontsize',14);
      set(gca,'fontsize',14);
     
      F_hhw_phi_ST_rea_mean = [age,Cons_HHW_rea_mean,Cons_PHI_rea_mean,Cons_ST_rea_mean];
      writetable(array2table(F_hhw_phi_ST_rea_mean),'./results/F_hhw_phi_ST_rea_mean_HHW_PHI_ST.xlsx')
            
     %%     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %%%%% CONSUMPTION - SD OF FIRST DIFFERENCES %%%%%%%%
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
      Cons_ST_abi =Wio_abi-Ppaid_spot;
      Cons_HHW_abi =Consumption_HHW(Cg_abi,iX);
      Cons_PHI_abi =Wio_abi-Ppaid_N;
      
      Cons_ST_rea =Wio_rea-Ppaid_spot;
      Cons_HHW_rea =Consumption_HHW(Cg_rea,iX);
      Cons_PHI_rea =Wio_rea-Ppaid_N;

      %%%%%%%%%%%%%
      % abi
      %%%%%%%%%%%%%
            
      Delta_Cons_ST_abi  = Cons_ST_abi(2:T,:)-Cons_ST_abi(1:(T-1),:);
      Delta_Cons_ST_abi_mean = sum(Delta_Cons_ST_abi.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
      Delta_Cons_ST_abi_sd=sqrt(sum(((Delta_Cons_ST_abi-repmat(Delta_Cons_ST_abi_mean,1,nInd)).*(iX(2:T,:)~=Nstates+1)).^2,2)./sum(iX(2:T,:)~=(Nstates+1),2));

      Delta_Cons_HHW_abi  = Cons_HHW_abi(2:T,:)-Cons_HHW_abi(1:(T-1),:);
      Delta_Cons_HHW_abi_mean = sum(Delta_Cons_HHW_abi.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
      Delta_Cons_HHW_abi_sd=sqrt(sum(((Delta_Cons_HHW_abi-repmat(Delta_Cons_HHW_abi_mean,1,nInd)).*(iX(2:T,:)~=Nstates+1)).^2,2)./sum(iX(2:T,:)~=(Nstates+1),2));

      Delta_Cons_PHI_abi  = Cons_PHI_abi(2:T,:)-Cons_PHI_abi(1:(T-1),:);
      Delta_Cons_PHI_abi_mean = sum(Delta_Cons_PHI_abi.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
      Delta_Cons_PHI_abi_sd=sqrt(sum(((Delta_Cons_PHI_abi-repmat(Delta_Cons_PHI_abi_mean,1,nInd)).*(iX(2:T,:)~=Nstates+1)).^2,2)./sum(iX(2:T,:)~=(Nstates+1),2));
     
   
      plot(age(2:T,:),Delta_Cons_HHW_abi_sd,age(2:T,:),Delta_Cons_PHI_abi_sd,age(2:T,:),Delta_Cons_ST_abi_sd)

      F_hhw_phi_ST_abi_Delta_SD = [age(2:T,:),Delta_Cons_HHW_abi_sd,Delta_Cons_PHI_abi_sd,Delta_Cons_ST_abi_sd];
      writetable(array2table(F_hhw_phi_ST_abi_Delta_SD),'./results/F_hhw_phi_ST_abi_Delta_SD_HHW_PHI_ST.xlsx')
      
      %%%%%%%%%%%%%
      % rea
      %%%%%%%%%%%%%
      
      
      Delta_Cons_ST_rea  = Cons_ST_rea(2:T,:)-Cons_ST_rea(1:(T-1),:);
      Delta_Cons_ST_rea_mean = sum(Delta_Cons_ST_rea.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
      Delta_Cons_ST_rea_sd=sqrt(sum(((Delta_Cons_ST_rea-repmat(Delta_Cons_ST_rea_mean,1,nInd)).*(iX(2:T,:)~=Nstates+1)).^2,2)./sum(iX(2:T,:)~=(Nstates+1),2));
      plot(age(2:T,:),Delta_Cons_ST_rea_sd)

      Delta_Cons_HHW_rea  = Cons_HHW_rea(2:T,:)-Cons_HHW_rea(1:(T-1),:);
      Delta_Cons_HHW_rea_mean = sum(Delta_Cons_HHW_rea.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
      Delta_Cons_HHW_rea_sd=sqrt(sum(((Delta_Cons_HHW_rea-repmat(Delta_Cons_HHW_rea_mean,1,nInd)).*(iX(2:T,:)~=Nstates+1)).^2,2)./sum(iX(2:T,:)~=(Nstates+1),2));
      plot(age(2:T,:),Delta_Cons_HHW_rea_sd)

      Delta_Cons_PHI_rea  = Cons_PHI_rea(2:T,:)-Cons_PHI_rea(1:(T-1),:);
      Delta_Cons_PHI_rea_mean = sum(Delta_Cons_PHI_rea.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
      Delta_Cons_PHI_rea_sd=sqrt(sum(((Delta_Cons_PHI_rea-repmat(Delta_Cons_PHI_rea_mean,1,nInd)).*(iX(2:T,:)~=Nstates+1)).^2,2)./sum(iX(2:T,:)~=(Nstates+1),2));
      plot(age(2:T,:),Delta_Cons_PHI_rea_sd)
                  
      plot(age(2:T,:),Delta_Cons_HHW_rea_sd,age(2:T,:),Delta_Cons_PHI_rea_sd,age(2:T,:),Delta_Cons_ST_rea_sd)

      F_hhw_phi_ST_rea_Delta_SD = [age(2:T,:),Delta_Cons_HHW_rea_sd,Delta_Cons_PHI_rea_sd,Delta_Cons_ST_rea_sd];
      writetable(array2table(F_hhw_phi_ST_rea_Delta_SD),'./results/F_hhw_phi_ST_rea_Delta_SD_HHW_PHI_ST.xlsx')

    %%  
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TABLE WITH CONTRACT TERMS    
    % first year equilibroum contract terms
    % premium and frontloading
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
      % ABI- HHW

    consumption = Cg_abi(1:Nstates,1)';
    first_year_costs = MU(1:Nstates,1)';
    premium = Wio_abi(1)-consumption;
    frontloading = premium-first_year_costs;    
    table11_hhw_abi = [premium;frontloading];
    
     % REA- HHW

    consumption = Cg_rea(1:Nstates,1)';
    premium = Wio_rea(1)-consumption;
    frontloading = premium-first_year_costs;
    table11_hhw_rea = [premium;frontloading];
    
    % ABI and REA - GLTHI
    premium = Pspot_N(1:Nstates,1)';
    frontloading = premium-first_year_costs;
    table11_ger = [premium;frontloading];

    table_terms = [first_year_costs;table11_ger;table11_hhw_rea;table11_hhw_abi];
    table_terms = round(table_terms,3);
    
    writetable(array2table(table_terms),'../../figsandtabs/table_terms.xlsx');

    %%
    %%%%%%%%%%%%%%%%%%%%%%%
    %%% LAPSATION %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%
 
      Delta_P_PHI = Ppaid_N(2:T,:)-Ppaid_N(1:(T-1),:);
      lapse_PHI = Delta_P_PHI~=0;
      lapse_PHI_mean = sum(lapse_PHI.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
      
      Delta_P_P_spot = Ppaid_spot(2:T,:)-Ppaid_spot(1:(T-1),:);
      lapse_Spot = Delta_P_P_spot~=0;
      lapse_PSPot_mean = sum(lapse_Spot.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
     
      lapse_HHW_abi = abs(Delta_Cons_HHW_abi)>0.000001;
      lapse_HHW_abi_mean = sum(lapse_HHW_abi.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 

      lapse_HHW_rea = abs(Delta_Cons_HHW_rea)>0.000001;
      lapse_HHW_rea_mean = sum(lapse_HHW_rea.*(iX(2:T,:)~=(Nstates+1)),2)./sum(iX(2:T,:)~=(Nstates+1),2); 
     
      F_lapse_rates =[age(2:T,:),lapse_PHI_mean,age(2:T,:),lapse_HHW_abi_mean,age(2:T,:),lapse_HHW_rea_mean,age(2:T,:)];
      writetable(array2table(F_lapse_rates),'./results/F_lapse_rates_PHI_HHWabi_HHW_rea.xlsx')
      
